arXiv:l502.07743V 1 [cs.SY] 21 Jan 2015 


Tracking an Object with Unknown Accelerations 
nsing a Shadowing Filter 

Kevin Judd, University of Western Australia 
March 2, 2015 


Abstract 

A commonly encountered problem is the tracking of a physical ob¬ 
ject, like a maneuvering ship, aircraft, land vehicle, spacecraft or animate 
creature carrying a wireless device. The sensor data is often limited and 
inaccurate observations of range or bearing. This problem is more difficult 
than tracking a ballistic trajectory, because an operative affects unknown 
and arbitrarily changing accelerations. Although stochastic methods of fil¬ 
tering or state estimation (Kalman filters and particle filters) are widely 
used, out of vogue variational methods are more appropriate in this track¬ 
ing context, because the objects do not typically display any significant 
random motions at the length and time scales of interest. This leads us 
to propose a rather elegant approach based on a shadowing filter. The 
resulting filter is efficient (reduces to the solution of linear equations) and 
robust (uneffected by missing data and singular correlations that would 
cause catastrophic failure of Bayesian filters.) The tracking is so robust, 
that in some common situations it actually performs better by ignoring 
error correlations that are so vital to Kalman filters. 


1 Introduction 

The literature on tracking of isolated or multiple objects in uncluttered and 
cluttered environments is broad and varied, with methodological approaches 
ranging from variational techniques to an extensive variety of Kalman filters, 
particles filters, and other sequential Bayesian filters [311171 HU]. 

Optimal tracking (filtering or state estimation) requires identification of the 
appropriate length and time scales of a object’s motion, and the magnitudes 
of the uncertainty (observational and dynamical noise) and the nonlinearity at 
these scales. There are three broad scenarios: (1) when the update cycle is 
fast, such in guidance control, Kalman filters are appropriate, because system 
behaviours are almost linear; (2) when nonlinearity is significant, and the noise 
largely dynamical, as opposed to observational, then particle filters are appropri¬ 
ate; and (3) when nonlinearity is significant, but the noise largely observational, 
then variational methods are most appropriate. Mathematically situations (1) 
and (2) best assume the underlying process is stochastic, whereas situation (3) 
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it is best to assume a deterministic dynamical system: the success and universal 
acceptance of stochastic methods, over the past decades, has lead to them being 
applied in situation (3), where they are not the most appropriate choice [5]. 

Here we present an approach to tracking in situation (3) based on shadow¬ 
ing filters, which derive from the modern concept of shadowing in dynamical 
systems theory; literally meaning to find a trajectory that shadows the obser¬ 
vations [6] . The methodology has its roots in the work of Laplace and Gauss 
fitting celestial orbits as curves [3], and subsequent least squares approaches 
for ballistic trajectories. Shadowing filters lie within the domain of variational 
methods, as in optimal control, but are subtly different [5]. Shadowing filters 
are not equivalent to 4D-variational assimilation often employed in meteolog- 
ical and oceanographic modelling and forecasting. Nor are they equivalent to 
dynamic programming approaches that finds a Viterbi path [ 5 . 

To develop our methodology, a one-dimensional, or scalar case, is considered 
first, which is then extended to the multi-dimensional vector case, where the 
observations are of the components of the Cartesian position vector. From this 
basis other relevant observation situations that are often encountered are consid¬ 
ered, for example, using range or bearing observations from one or more sensors. 
A significant problem that arises in these types of observation networks is the 
way the covariance of observational errors varies with position, in particular, the 
singularities that occur when the target and sensors are co-linear or co-planar. 
Possibly the most surprising outcome of the shadowing filter approach is that 
singular covariances of observations are not a significant problem, and targets 
can be tracked through missing data and singularities relatively easily; indeed, 
sometimes covariances can be ignored entirely, obtaining very efficient tracking 
filters. 

To keep the exposition of the algorithm and its benefits clear, we restrict 
attention to tracking an isolated vehicle in an uncluttered environment. It 
should be clear, once the methodology is understood, that since the shadowing 
filter assumes a tracked object maintains a contiguous trajectory it will perform 
well with multiple targets in cluttered environments. 


2 Formulation and implementation 

2.1 Scalar case 

Our initial interest is tracking the position of a point object in one dimension, 
given a sequence of noisy observations. Let £ K be the observed position at 
time ti for i = 0,... ,n, and af be the variance of the observational error. The 
object’s dynamics are modelled by its position £ R and velocity Vi £ R at tj, 
and constant acceleration £ M for ti < t < ti+i. For notational convenience, 
define n = ti+i — f. Our goal is to have pi close to Vi subject to the accel¬ 
erations not being excessively large or changeable. We might therefore choose 
to minimise the total square error ~ Pi)^ subject to accelerations 

over the interval being bounded, of < i = 0,..., n — 1. Bounded acceleration 
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is an appropriate constraint, but it introduces technical difficulties. Although 
these difficulties can be over-come [ 8 ], it is more convenient and efficient to in¬ 
stead constrain the root mean squared acceleration over the entire trajectory, 

Assuming Newton’s laws and Galilean transforms apply to the point object’s 
motion, then the stated optimisation problem can be posed using a Lagrangian: 


n 

n— 1 ^ 

+ X! -Pi- 

n—1 

+ ^ - Vi- atn) 

i =0 




- (4 



( 1 ) 

( 2 ) 

(3) 

(4) 


where Ai G K, G R and 77 > 0 are dual variables. Solving the optimisation, 
and defining 


P[t) =Pr+ Vi{t 


ti) + 2“i(^ “ 


for ti<t<U+i, 


(5) 


provides an optimal quadratic spline estimate of a particle’s path assuming 
piecewise constant accelerations. The optimal solution occurs where all the 
partial derivatives of L are zero: 
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Equation ((O]) is defined for 0 < i < n while (ITHTOI) are defined for 0 < i < n — 1. 
With the exception of (HU, the remaining five equations are linear in the un¬ 
knowns. However, ^ and ry are related through term (|3|) of the Lagrangian, which 
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is the only place they appear. Hence, one can solve the linear equations (EHini) 
for a fixed ij, then compute the corresponding value of ^ from eq. dTTI) . The op¬ 
timal solution for any ^ can be approximated arbitrarily closely by an efficient 
one-dimensional search, such as Brent’s method m- In practice it is unlikely 
that ^ needs to be specified precisely, after all, it is only a bound on the root 
mean squared acceleration. Consequently, it is usually sufficient to work only 
with 77 , treating it as a smoothing, or regularisation ^parameter. 

Combining m and m to eliminate the Vi giveqj 

Pi+iTi-i -piin + +pi-iTi = ^{uiTi + ai_iTi_i)n_iTi, ( 12 ) 

for 0 < z < n. Combining ®, © and © and eliminating the dual variables 
Xi and Pi (as explained in the following), obtains another set of expressions 
relating the Pi and , which when combined with ( 1121 ) enable solving for a near 
optimal solution very efficiently. 

There is a certain amount of redunancy in the equations just stated, but 
to assist formulation of a solution using matrix notation it is advantageous to 
retain the redundancy. 

Define column vectors V = (Vo, ■ ■ ■ p = {po,... ,pn)^ £ 

]R"+\ A = (Ao,..., A„_i)^ G R”, p = {po ,..., Pn-iV G 1^": and finally a = 
(oo,..., a„_i)^ G R.". Define r to be the nxn matrix of zeros with main 
diagonal (tq, ..., Tn-i), and I to be the (n -I- 1) x (n -I- 1) matrix of zeros with 
main diagonal ,..., Define a n x n matrix D to have all entries zero 

except — 1 on the main diagonal and 1 on the first lower diagonal, and similarly 
define a (n -I- 1) x n matrix E\ specifically 

{ - 1 , *=j, 

1 , *=j + l, (13) 

0 , otherwise, 

when the entry is defined. It follows that the linear equations ®, © and ® 
can be succinctly expressed as 


EX=I{V—p), Dp = tX, 


2riTa 


1 


ttX + Tp. 


(14) 


In the last of the three equations of m, T is invertible, so that a r factor can 
be canceled on the left of each term. 

Define a nxn matrix L, and nx (n-|-1) matrix M, to have a lower triangular 
form: 


Lij — ^ij — 


- 1 , 

0 , 


i > j, 

otherwise. 


(15) 


^To do this multiply m by Ti—i, then take another copy of m with i replaced with i — 1, 
multiply by r^, and substract this from the former. Then use (IIOII to eliminate Vi — Vi—i. 
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when the entry is defined. It can be easily verified that DL = I and EM = J 
where I is the n x n identity matrix, and J is a (n + 1) x (n + 1) matrix with 


i l, i = j and i ^ n + I, 

— 1, i = n + 1 and j n + 1, (16) 

0, otherwise. 

The identities DL = I and Dfi = tX impljH fj, = LtX. The identities EM = J 
and EX = I{V — p) impljH that A = MI{V — p) and — Pi) = ti¬ 

lt follows, by substitution into the last equation of (ITTl) . that 

2r]a = {^T-M + LrM^ X{V - p), (17) 

which can be combined with (|T^ as follows. Define (n — 1) x (u + 1) matrices 
A and B, and (n — 1) x n matrix G: 


{ rfn+i, 

0 , 


i=j, 

i + l=j, 
otherwise, 


(18) 


'^ 2+1 5 - jf 

-{Ti + n+i), 1 + 1 = j, 

n, 1 + 2 = j, 


0, otherwise. 


(19) 


A=^G (^^T-M + Ltm'^ . ( 20 ) 

Then (TT^ can be written Ep = ^Ga, which when combined with (ITTl) obtains 
the equation 

(AX + r]B)p = AXV, (21) 

and the additional constraint YYo ~ Pi) = 0. For stability reasons 

discussed in section 1^31 this constraint will be extended to Yl)=o = 

0. Defining a n x (n +1) matrix S to be i? augmented with a final row of zeros, 
and a n X (n + 1) matrix A to be A augmented with a final row of ones, then 


(AX + riB)p = AXV, (22) 

encodes both m and the extended constraint. Solving ([2^ by singular value 
decomposition obtains a least squares approximate solution for the position p 
for given smoothing parameter rj. See section [2.31 on the nature of this approxi¬ 
mation and the optimal presentation of the data in V. It is important to read 
section\2^ before implementing ISSI) . 

^If DL = I, then DLt\ = tA, but Dfj, = tA, implying n = LtX. 

®If EM = J, then EMXIV — p) = JX{P — p). If EX = XfP — p) also, then, since 
the first n rows of J is the n X n identity, A = MX('P — p). Substituting this A back into 
EX = X{V — p), gives JX(P — p) = XI(P — p); the first n rows are a tautology, but the last 
implies YiEo - Pi) = 0. 
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2.2 Vector case 


Consider now the situation where a point vehicle is positioned in a d-dimensional 
Euclidean space with Cartesian coordinates. Suppose that the coordinate posi¬ 
tions are observed as a sequence Vi S in such a way that the observational 
errors have a d x d covariance matrix Ci, and corresponding information matrix 
li = C~^. The quantities to be determined are pi G , Vi G , and at G 
which are all now d-dimensional column vectors. If the aim is to track the tra¬ 
jectory under the assumption of bounded RMS magnitude of the acceleration, 
the Lagrangian o now becomes vectorised as 

1 " 

L = -;^^{'Pi-Pr)'^Ti[Vi-pr) 

n—1 ^ 

+ ^ - Pi- V^n - 

i^O 
n—1 

-b ^ pf+i{vi+i - Vi- a^n) 
i^O 

n-1 \ 

^ ^ "^i^i i^n I : 

i=0 ) 

where Xi G , pi G R^^ are now d-dimensional column vectors, but 77 G R. The 
superscript T indicates the transpose. 

The solution of the vectorised optimization problem proceeds identically to 
the scalar case, using the same linear algebra methods. Let V G denote 

a column vector being the time-series of n -I- 1 observations Vi G R'^ stacked in 
d-dimensional blocks, and similarly for position variables p G Let I 

denote the (n + l)d x (n -I- l)d block diagonal matrix with the d x d information 
matrices Zi along the diagonal. Finally, let Id denotes the d x d identity matrix 
and let M = M® Id denote the outer product of an arbitrary matrix M with Id, 
that is, M has a block structure where each scalar entry of M becomes a 
d X d-block Mijid of M. Then the vectorised solution is 

(AI + rjBj p = MV. (27) 

2.3 Implementation and interpretation of the filter 

Solving or ((?7l) obtains an approximate solution to the optimal shadowing 
trajectory for a given smoothness rj. To see this, note that the system of equa¬ 
tions (I22[) is under-determined: there are n linear equations in n -|- 1 unknowns. 
This occurs because in deriving (I12|) the velocity variables were eliminated, but 
to completely define a trajectory the velocity needs to be known at some time; 
usually the initial or final velocity is specified or solved for. To solve for the 
velocity requires introducing another n variables and n equations to solve for 



(23) 

(24) 

(25) 

(26) 
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all the velocities, which is significant additional computation for very little ben¬ 
efit. The approximation (1221) relies on the fact that if the time window of the 
trajectory is sufficiently long, then accurate specification of the initial velocity 
is not required. It just means the initial part of the trajectory may not accu¬ 
rately fit the observations. However, leaving the initial velocity unspecified can 
lead to instability for short time windows. Imposing the extended constraint 
Sr=o ~ Pi) = 0 overcomes possible instability by implicitly defining an 

initial velocity. 

In formulating the Lagrangian forward differences were used to express po¬ 
sition in terms of velocity and acceleration. Unfortunately, this results in A 
and B having a lower triangular form. Consequently, the approximation er¬ 
rors are largest for the pi with largest z, which is not what is wanted for state 
estimation and forecasting; it is preferable that the smallest errors are at the 
most recent times. Reformulating the Langrangian with backward differences 
solves this problem, however, there is much simpler solution: initially reversing 
the time-series data sequence (ri,Pi,Ii), applying the filter (1^ or (l?71) . then 
reversing shadowing trajectory time-series pi to obtain the desired result. Even 
this trick is unnecessary. Let R denote the matrix that reverses a vector, then 
the time-series reversal trick, is equivalent to changing m to 

{AI + pB) Rp = ARXV, (28) 

but since R~^ = R, multipling on the left by R obtains 

{RARJ + pRBR) p = RARIV, (29) 

where the matrices RAR and RBR are just A and B with their rows and columns 
reversed. Hence, the time-reversal trick is some bookkeeping when constructing 
A and B. 


3 Illustrative examples 

This section provides demonstrations of the use of the proposed methods. The 
scalar filter is considered first, both as an off-line smoothing filter and sequential 
state-estimator. The vector filter is considered for observations in cartesian 
coordinates, with and without correlation, which is a preliminary to section |T] 
where non-cartesian observations are considered. 

3.1 Scalar filter for moothing and sequential tracking 

Here tracking of one observed variable is examined for increasing values of 
smoothing paramter p. The filter is employed as smoothing filter over the en¬ 
tire observation window, and as a sequential state-estimator. In both cases the 
time-reversal trick discussed in section 12.31 is employed. In this example em¬ 
ploys a large red noise component Xt to mimic a vehicle maneuvering in an 
unpredictable way. 
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Time 


Figure 1: Position tracking 25 +10 sin(t/15) +xt + 3et, 0 < t < 100 where et is a 
white noise process A^(0,1), and Xt an independent cumulative of a white noise 
process A^(0,1). Shadowing filter results for smoothing parameters as stated. 



Figure 2: Computed accelerations for tracking shown in fig. [T] for selected f]. 
Accelerations are scaled by y/rj to allow easier comparision. 
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Time 


Figure 3: Sequential tracking of same data as shown in fig. [T] for selected rj. 
State estimates use only observations up to that time, and final filtered state at 
that time is plotted. 
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Figured] reveals how the implied approximation in the solution (1^^ results 
in the position tracking deviating from ideal at the beginning (t = 0 ) of the 
time series. Without the time-reversal trick, this deviation would have occurred 
at the end (t = n); although the deviation is small, it is significant, but if the 
time-series window is large enough, then there is no significant effect for t > 10 . 
Optimal smoothing appears occur in the range 10 < ry < 100. 

Figure [2] shows how smaller ry result in large and rapidly switching acceler¬ 
ations, while larger rj result in much smaller accelerations applied over longer 
periods. 

Figure El demonstrates using the same filter as in figure dl as a sequential 
state-estimator; the filter is applied only to the observations up to that time. In 
the this tracking mode the smoothing parameter 77 is seen to act like an inertial 
damping. For the larger rj = 10000 the tracking lags the true trajectory. For 
smaller rj values the tracking is better, but note how a sequence of observations 
with repeated negative bias for 60 < t < 75 result in the tracking over-shotting 
the turn near t = 70. When the repeated bias ends around t = 75 the near 
optimal 77 = 100 track jumps back to good estimates, while the 77 = 1000 track 
turns back smoothly toward the true trajectory. 

3.2 Vector filter with uncorrelated observations 

If the observations of each component of the d-dimensional position are uncor¬ 
related, then filtering can be accomplished very efficiently using a scale filter, 
which will come in useful later when non-cartesian observations are considered. 

When the observations of each component of the d-dimensional position are 
uncorrelated the covariance matrices Ci are all diagonal, and it is unnecessary 
to use the vectorised hlter (l27)) . which has been expanded by an outer product 
with Id; instead it is sufficient to solve (l2^ seperately for each component. If 
all the components have proportionally the same variance at each time, then the 
singular value decomposition only needs to be solved once for all components, 
that is, one information matrix I is needed, whose elements are proportional 
to the variances, and V becomes a (n -|- 1) x d matrix of observations, so that 
eq. (EH) then solves for all components of p simultaneously. 

Figure m shows tracking in two dimensions in this situation. Observe how 
a very large 77 = 10000 results in poor tracking as the tracking curve is pulled 
toward the mean of the observations. The poor tracking at the beginning can 
also be observed for larger 77 . Fig|4jb) shows the important case of sequential 
estimates, that is, sequential tracking of the object using all the observations 
obtained up to a given time. For efficiently reasons, one would in practice only 
use a finite window of past observations. Details of how to determine the optimal 
window for a given system and purpose is beyond the scope of this paper and 
is discussed in the general context of shadowing filters elsewhere Ha¬ 


lo 


(a) (b) 




Figure 4: Position tracking for 0 < t < 150 of the path {x, y) = 10(t — 10)/150 + 
(1/3)(1 — t)(sin(t/15), 2 —1/15). Observational errors independent on each com¬ 
ponent with standard deviation 5. (a) Final tracking curves for various 77 . 

(b) Sequential position estimates for 77 = 1000, that is, last state of a shad¬ 
owing trajectory obtained using all observations up to a given time. 
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4 Using non-Cartesian observations 

A frequently encountered tracking problem involves using bearing and range ob¬ 
servations, or combinations of multiple bearing or range observations. Several 
important situations are worth considering. Active radar location uses range 
and bearing information. Satellite interferometry uses range and bearing in¬ 
formation, but the range is much more accurately measured than the bearing. 
Global positioning by satellite uses only range information, but from multiple 
reference satellites. Tracking wireless devices can use range information inferred 
from signal power at multiple transponders. Passive sonar location provides 
bearing information, but poor range; often bearings from multiple sensor points 
are used. 

All of the applications mentioned can be dealt with using the vector filter 
(Ell), by transforming the observations into raw Cartesian position estimates and 
computing the appropriate information matrix. The transformations are simple 
geometry, but computing the information matrices requires some approximation 
or restrictions. 

Let p G he the position in Cartesian coordinates, and let q G he 
a vector of d noise-free observations of the position is some other coordinates. 
Suppose there is an invertible function /, on some domain, such that p = f{q). 
Given a noisy observation Q — q + Ag, the transform / provides a raw position 
estimate V = f{Q) = p -\- Ap. Given the covariance matrix Cq of the observa¬ 
tions, the covariance Cp of the estimate is required. The column vector Ag is the 
error in the observation, and to a first approximation, the error in the estimate 
is Ap « TAg, where J = dqf{q) is the Jacobian matrix of / at g. Since the 
covariance Cp is the expected value of the outer product ApAp^, it follows that, 
to a hrst approximation, Cp = JCqJ^. It also follows that the corresponding 
information matrices are related by Ip = K^IqK, where K = J~^ = dp{f~^). 
Note that since only the information matrix is needed in the shadowing filter 
under discussion, it is sometimes easier to compute the K directly using f~^, 
than it is to compute J and invert it. This is the case in some of the following 
examples. 

An important problem, which will be returned to in each of the following 
sections, is that although the correlation matrix Cq may be well known, the 
transformation matrices J and K depend on the target’s location, which is un¬ 
known. The raw position estimate V = f{Q) could be used, but this introduces 
errors in the supposed covariance. When the transformed coordinates are highly 
correlated and the transform very non-linear, then a small error in the raw po¬ 
sition estimate can give rise to a very wrong estimate in the correlation. This 
problem plagues Kalman filters, and other filters that need a covariance esti¬ 
mate to reliably estimate the state. A signihcant advantage of a shadowing filter 
is the robustness gained from finding a shadowing trajectory, rather than just 
a current position estimate. This robustness means that the correlation in the 
raw position estimates are of little importance, that is, ignoring the correlation 
can have little effect on the quality of the tracking. 
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4.1 Range and bearing observations 

Consider a target tracked in the plane, position p = {x,y), using observations 

q = {r,9), where r is the range from a reference point (a, 6) and 0 the bear¬ 

ing in radians measured in the anti-clockwise direction from the x-axis. The 
transformation p = f{q) is given by 

X = a + r cos 6, (30) 

y = b + rsinO. (31) 

Under the assumption that r is not close to zero, and the variances of r and 9 are 
small, then the covariance and information matrix of x and y are approximated 
as previously described using 

r_fcos9 -rsm9'\ , . 

\^sin0 rcos9 j 
or 

cos 6* sin6» \ . 

\^—(1/r) sin6* (l/r)cos0y’ 

Figure [5] shows the tracking of a target using range and bearing information 
of different accuracy. In panels (a) and (b) the correlation of the raw position 
estimates is ignored, and the results are good, that is, the shadowing filter pro¬ 
vides significant improvement over the raw position estimates. Panel (c) uses the 
same data as panel (b), but tries to take into account the correlation of the raw 
position estimates using the correlation computed at the raw position estimate; 
the result is worse than assuming no correlation. If the same is attempted for 
the data of panel (a), the result is a worse failure, because the radius is poorly 
estimated and since this appears as a reciprocal in the transformation, small 
errors in the radius can lead to very poorly estimated correlations, so much so 
that the quality of the filtering is much worse. 

This example provides an excellent illustration of how the robustness of a 
shadowing filter has significant gains over other filters. Not only does ignoring 
the correlation result is better tracking, it is also more efficient, because rather 
than using the vectorised filter (EZl), the simpler, more compact, scalar filter (El]) 
can be used. 


4.2 Multiple bearing observations 

Consider a target tracked in the plane, position p = (x, y), using bearing obser¬ 
vations q = {9,9') from two distinct reference points (a, 6) and {a',b'). Under 
most circumstances there exists unique s, s' G R such that 

{x,y) = {a,b) + s{cos9,sm9) (34) 

= {a ,b') + s'{cos9',sm9'). (35) 


Hence, a raw position estimate can be found by solving the linear equations 


/cos 9 — cos 9'\ f s\ _ fa' — a\ 

l^sin0 -sin 07 \^s7 “ Vfe'-fej ’ 


(36) 
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Figure 5: Position tracking of same path as fig. ID Circle with cross-hair is 
the observation site. Sequential position estimates and final trajectory estimate 
for T] = 1000. (a) Bearing measurement ten times more accurate than range, 
(b) Range measurement ten times more accurate than bearing. In both (a) and 
(b) the estimation ignores correlations and treats each component as being in¬ 
dependent. (c) As case (b) but using correlation as computed from raw position 
estimate; this does worse than (b), because the correlation is calculated relative 
to the raw position estimates, which are very misleading. 
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provided the equations are consistent and non-singular. Since the angles are 
observed angles, inconsistency can occur when d' ~ ±0. For some sequential fil¬ 
ters such nearly-singular situations can be devastating, but, since the shadowing 
filter is estimating a trajectory from a sequence of observations, there is gener¬ 
ally no harm in simply dropping observations corrupted by near-singularities, 
or replacing them with forecasted positions, or crudely interpolated positions. 
There is generally no harm in doing either of these for short periods. Dropping 
observations requires using a larger r time-gap between observations. If the 
observations were equally spaced in time, this leads to a lot of special compu¬ 
tation, and so in this case it is generally easier to insert a forecasted position 
with a suitably scaled down information matrix to account for the errors in the 
forecast; see details given later. 

In this multiple-bearing situation it is difficult to compute J directly, but K 
is easily computed. Note that 


Q — arctan 


y-h 


and 6’ = arctan 


X — a 


y-b' 


X — a' 


(37) 


It follows that 


({x-a)/r'^ {y-b)/r^\ 

{ix-a')/r'^ {y-b')/r'y’ 


(38) 


where r'^ = {x — a)^ + (y — b)^ and = {x — a'Y + (y ~ b')"^^ which requires 
only trivial computation once an estimate of (cc, y) is obtained. 

As mentioned, a significant problem arises when the sensors and target are 
collinear, 9 k. ±0', because the linear equations (1361) are singular or badly condi¬ 
tioned. This can result in raw position estimates far their true position. Figure[6] 
shows an example of tracking using two mobile sensors for detecting bearing, 
and a mobile target. In this example the target moves on a circular path clock¬ 
wise from the 12 o’clock position. When the target is between the 4 and 5 
o’clock position it is directly between the sensors, and at the end of its path, 
around the 7 o’clock position, the target is almost directly behind both sensors. 
Both of these situations lead to a poorly conditioned matrix in eq. (1361) and 
hence poor raw position estimates. 

For the estimates shown in fig. [5] the component correlations of the raw 
position estimate are ignored, as in fig [5Da) and (b). The ill-conditioning is 
dealt with by using the 1-norm estimate of the reciprocal condition as returned 
by LAPACK [2]. This number varies between zero and one, with small values 
indicating bad conditioning. This norm is a natural candidate to scale the 
information matrices to give raw position estimates from poorly conditioned 
situations small weight. 

From fig. [6] it can be seen that when the raw position estimates are obtained 
under good conditioning, the sequential estimates (large circles) are good. Under 
poor conditioning (small circles) the sequential estimates are mainly forecasts 
from the preceding trajectory positions, and the wild raw position estimates are 
ignored. When the target passes beyond the 5 o’clock position and conditioning 
improves, and the sequential estimates return to good position estimates. Note 
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Figure 6: Position tracking using bearings from two moving sensors. For 0 < f < 
100 one sensor moves between (—3, 3) and (3,1), and the other between (—3, —2) 
and (3, —1), both at constant speed, while the target moves on the circular path 
(sin(t/25), cos(t/25)). The centre of the circles is the position of the sequential 
state estimates, the diameter of the circle indicates the condition-number weight 
applied to that position’s raw position estimate. 
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also that the final trajectory almost exactly matches the true trajectory through 
the 4 to 5 o’clock position. 

4.3 Multiple range observations 

Consider a target tracked in the plane, position p = {x,y), using range obser¬ 
vations q = (r,r') from two distinct reference points (a, 6) and {a',b'). Un¬ 
der most circumstances the location can be obtained from the solutions of 
r'^ = [x — aY + {y — bY and r'^ = (x — a'Y + {y — b'Y, assuming that the 
non-uniqueness can be resolved. By taking partial derivatives of these two 
equations with respect to x and y, it follows implicitly that 



(39) 


5 Partial observations 

Our stated formulation of the tracking problem allows for non-uniformly spaced 
observations, but all examples thus far have used only uniformly spaced obser¬ 
vations. We make a simple demonstration using non-uniformly spaced observa¬ 
tions by considering a situation where observations are missing. 

Figure [7] demonstrates tracking using similar data to figure [U where 75% of 
the observations are missing, comparing this to the tracking calculations when 
all observations are available. Two values of the smoothing parameter r/ are 
used. These results demonstrate that the tracking algorithm is very robust. 
The position tracking is very similar when there is missing observations. The 
acceleration estimates are also very similar. Overall the tracking is slightly 
smoother, and accelerations less variable, when there is missing data, but this is 
something of an artifact, because the effective amount of smoothing for a fixed 
rj depends on the amount of observations available; less data results in more 
smoothing. 

6 Conclusion 

Under the assumption of piecewise constant accelerations a shadowing filter 
algorithm has been derived and implemented efficiently for scalar time-series 
of observations. Vector time-series of observations can be dealt with efficiently 
using the same algorithm if each component is observed with uncorrelated errors. 

The scalar algorithm can be easily extended to deal with tracking in d- 
dimensions where the position vector components are not observed directly and 
a transformation of the observations can be used to obtain an initial raw position 
estimate. The question then arises of how to deal with the correlation of the 
errors this introduces. Remarkably, experiments reveal, as illustrated in figs 
H M and 1 that ignoring this correlation has little significant effect. Simply 
using the raw position estimates to estimate the correlations produced worse 
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Figure 7: Position tracking and accelerations using similar data to figure [T] and [2l 
but with partial observations for two rj smoothing values. The first two curves 
show calculations when 75% of observations are missing, with crosses marking 
the missing observations, and the second two curves using all observations. 
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position estimates, because errors in the raw position estimates give misleading 
indications about the correlation. This problem is unvoidable to Kalman filters. 
It may be that some more complex algorithm could be devised to better estimate 
the correlations, but this will increase the amount of computation, for possibly 
no significant gain. Just taking correlations into account in the stated algorithm 
increases the size of matrix requiring singular value decomposition from n(n+l), 
to (Pn{n + 1) entries, so the computation cost is significant. 

There are a number of other implementation issues that have not been dis¬ 
cussed, the most important of which is the optimal window size n for obtaining 
a shadowing trajectory and position estimates. The window size is problem de¬ 
pended, but a method for determining an appropriate window size is discussed 
at length elsewhere [12]. This cited work also discusses issues of how best to 
implement sequential filtering. 
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